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Abstract 

A physical-mathematical approach to anomalous diffusion may be based 
on generalized diffusion equations (containing derivatives of fractional order 
in space or/and time) and related random walk models. By the space- 
time fractional diffusion equation we mean an evolution equation obtained 
from the standard linear diffusion equation by replacing the second-order 
space derivative with a Riesz- Feller derivative of order a € (0, 2] and 
skewness < min{a,2 — a}), and the first-order time derivative with 
a Caputo derivative of order j3 G (0, 1] . The fundamental solution (for the 
Cauchy problem) of the fractional diffusion equation can be interpreted as 
a probability density evolving in time of a peculiar self-similar stochastic 
process. We view it as a generalized diffusion process that we call fractional 
diffusion process, and present an integral representation of the fundamental 
solution. A more general approach to anomalous diffusion is however known 
to be provided by the master equation for a continuous time random walk 
(CTRW). We show how this equation reduces to our fractional diffusion 
equation by a properly scaled passage to the limit of compressed waiting 
times and jump widths. Finally, we describe a method of simulation and 
display (via graphics) results of a few numerical case studies. 
Mathematics Subject Classification 2000: 26A33, 33E12, 33C60, 44A10, 
45K05, 47G30, 60G18, 60G50, 60G51, 60G55, 60J60, 60J70. 
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Physics, No. 621]. The present e-print is a revised version (with up-date annotations and 
references) of that contribution, but essentially represents our knowledge of that early 
time. 



1 Introduction 



It is well known that the fundamental solution (or Green function) for the 
Cauchy problem of the linear diffusion equation can be interpreted as a 
Gaussian (normal) probability density function (pdf) in space, evolving in 
time. All the moments of this pdf are finite; in particular, its variance 
is proportional to the first power of time, a noteworthy property of the 
standard diffusion that can be understood by means of an unbiased random 
walk model for the Brownian motion. 

In recent years a number of master equations have been proposed for 
random walk models that turn out to be beyond the classical Brownian 
motion, see e.g.Klafter et al. [34]. In particular, evolution equations 
containing fractional derivatives have gained revived interest in that they are 
expected to provide suitable mathematical models for describing phenomena 
of anomalous diffusion, strange kinetics^ and transport dynamics in complex 
systems. Recent references include e.g. [3]l2igi2HI23E9ESI39l32l231l3Bl 

ED EI]. 

The paper is divided as follows. In Section 2 we introduce our fractional 
diffusion equations providing the reader with the essential notions for the 
derivatives of fractional order (in space and in time) entering these equations. 
More precisely, we replace in the standard linear diffusion equation the 
second-order space derivative or /and the first-order time derivative by 
suitable integro- differential operators, which can be interpreted as a space or 
time derivative of fractional order a G (0, 2] or f3 G (0, 1] , respectively. The 
space fractional derivative is required to depend also on a real parameter 
9 (the skewness) subjected to the restriction \9\ < min{a, 2 — a}. Then, 
in Section 3 we pay attention to the fact that the fundamental solutions 
(or Green functions) of our diffusion equations of fractional order in space 
or/and in time can be interpreted as spatial probability densities evolving 
in time, related to certain self-similar stochastic process. We view these 
processes as generalized (or fractional) diffusion processes to be properly 
understood through suitable random walk models that have been treated by 
us in previous papers, see e.g. [151 IS CEH1 12Ql [211 E3]- In Section 4 we show 

2 To the topic of strange kinetics a special issue (nowadays in press) of Chemical Physics 
is devoted where the interested reader can find a number of applications of fractional 
diffusion equations 

3 We remind that the term "fractional" is a misnomer since the order can be a real 
number and thus is not restricted to be rational. The term is kept only for historical 
reasons, see e.g. [17]. Our fractional derivatives are required to coincide with the standard 
derivatives of integer order as soon as a = 2 (not as a = 1!) and j3 = 1 . 
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how such evolution equations of fractional order can be obtained from a more 
general master equation which governs the so-called continuous time random 
walk (CTRW) by a properly scaled passage to the limit of compressed 
waiting times and jump widths. The CTRW structure immediately offers a 
method of simulation that we roughly describe in Section 5 where we also 
display graphs of a few numerical case studies. Finally, in Section 6, we 
point out the main conclusions and outline the direction for future work. 

2 The space-time fractional diffusion equation 

By replacing in the standard diffusion equation 
d d 2 

—u(x,t) = ——^u(x,t), — oo < x < +oo , i>0, (2.1) 
at ox 1 

where u = u(x,t) is the (real) field variable, the second-order space 
derivative and the first-order time derivative by suitable integro- differential 
operators, which can be interpreted as a space and time derivative of 
fractional order we obtain a generalized diffusion equation which may be 
referred to as the space-time-fractional diffusion equation. We write this 
equation as 

t D^u(x,t) = x D%u(x,t), -oo<x<+oo, t>0, (2.2) 

where the a , 6 , [3 are real parameters restricted as follows 

0<a<2, \6\ < min{a,2 - a} , < (3 < 1 . (2.3) 

In (2.2) x Dq is the space-fractional Riesz-Feller derivative of order a and 
skewness 9 , and jD* is the time-fractional Caputo derivative of order j3 . 
The definitions of these fractional derivatives are more easily understood 
if given in terms of Fourier transform and Laplace transform, respectively. 
Generically, u(x, t) is interpreted as a mass density or a probability density 
depending on the space variable x, evolving in time t. 

In terms of the Fourier transform we have for the space-fractional Riesz- 
Feller derivative 

F{ x D$f(x);K} = -<(k) /(«), <(„) = l^e^ 11 ^/ 2 , (2.4) 

where f{n) = T {f(x); k} = e +i-nx ^ x _ j n ^]j er wor( j s the symbol 
of the pseudo-differential operator x Dq is required to be the logarithm of the 
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characteristic function of the generic stable (in the Levy sense) probability 
density, according to the Feller parameterization |12^ [13] . We note that the 
allowed region for the parameters a and 9 turns out to be a diamond in 
the plane {a, 9} with vertices in the points (0,0), (1,1), (2,0), (1,-1), 
that we call the Feller- Takayasu diamond^. For a = 2 (hence 9 = 0) we 



have T { X D% f(x)\K} 



{—in) , so we recover the standard second 



derivative. More generally for 9 = we have f { x Dg f(x); k} = —\n 

-( K 2 ) a / 2 SO 

' l dx 2 J 



(2.5) 



In this case we refer to the LHS of (2.5) as simply to the Riesz fractional 
derivative of order a . Assuming a / 1,2 and taking 9 in its range one can 
show that the explicit expression of the Riesz-Feller fractional derivative 
obtained from (2.4) is 



where, see [18], 

sin [(a - 9) vr/2] 



(a, 9) 



[c + (a,9) x Dl + c-(a,9) x D a _] f(x), 
sin [(a + 6)n/2] 



sin (an) 



c-(a,9) 



and X D± are Weyl fractional derivatives defined as 



*D%f(x) 



dx 

£_ 

dx 2 



jl—ct 



f(x) 



T 2-a 



f^) 



if 
if 



sin(a7r) 



< a < 1 



1< a < 2. 



(2.6) 



(2.7) 



(2. 



In (2.8) the X I± (/i > 0) denote the Weyl fractional integrals defined as 

Ox > 0) (2.9) 



m j- 
i 



(t-xr^KOda. 



oo 



4 Our notation for the stable distributions has been adapted from the original one by 
Feller. From 1998, see [18] , we have found it as the most convenient among the others 
available in the literature, see e.g. [32, 36, 45, 47, 53, 54, 61 . Furthermore, this notation 
has the advantage that the whole class of the strictly stable densities is represented. As 
far as we know, the diamond representation in the plane {a, 9} was formerly given by 
Takayasu in his 1990 book on Fractals 59 . 
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In the particular case 6 = we get c+(a,0) = c_(a,0) = l/[2 cos (an/2)] , 
and, by passing to the limit for a — > 2~ , we get c+(2, 0) = c_(2, 0) = —1/2 . 
For a = 1 we have 



x D$f(x) 



cos(^vr/2) + sin(#7r/2) X D f(x) 



(2.10) 



where f(x) = — f(x) , and 



x Dl Six) 



d_ 

dx 




•+oo m 



da- (2.11) 



— oo x 4 



In (2.11) the operator X H denotes the Hilbert transform and its singular 
integral is understood in the Cauchy principal value sense, see [20] . 

The operator x Dg has been referred to as the Riesz-Feller fractional 
derivative since both Marcel Riesz and William Feller contributed to its 



Let us now consider the time-fractional Caputo derivative. Following the 
original idea by Caputo [6], see also [5j [T71 09] , a proper time fractional 
derivative of order j3 G (0, 1) , useful for physical applications, may be defined 
in terms of the following rule for its Laplace transforrrH 



c{ t D^f(t);s} = s /3 f(s)-s /3 - 1 f(0 + ), 0<P<1, (2.12) 

where f(s) = £{f(t);s} = e~ s ^ f(t) dt . Then the Caputo fractional 
derivative of f(t) turns out to be defined as 



In other words the operator is required to generalize the well-known rule 
for the Laplace transform of the first derivative of a given (causal) function 
keeping the standard initial value of the function itself. 

The space-time fractional diffusion equation (2.2) contains as particular 
cases the strictly space fractional diffusion equation when < a < 2 and 

5 Originally, in the late 1940's, Riesz |50| introduced the pseudo-differential operator 
x Iq whose symbol is \n\~ a , well defined for any positive a with the exclusion of odd 
integer numbers, afterwards named the Riesz potential. The Riesz fractional derivative 
x Dq := — xIq°" defined by analytical continuation was generalized by Feller in his 1952 
genial paper [T5] to include the skewness parameter of the strictly stable densities. 

6 For our purposes we agree to take the Laplace parameter s real 

7 The reader should observe that the Caputo fractional derivative differs from the usual 
Riemann-Liouville fractional derivative which, defined as the left inverse of the Riemann- 



definitiorjl. 




< p < 1. 



(2.13) 
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(3 = 1, the strictly time fractional diffusion equation when a = 2 and < 
(3 < 1 , and the standard diffusion equation (2.1) when a = 2 and = 1. 
For the equation (2.2) we consider the Cauchy problem 

u(x, + )=ip(x), xGR, u(±oo, i) = , t>0, (2-14) 

where tp(x) is a sufficiently well-behaved function. By its solution we 
mean a function u a g(x, t) which satisfies the conditions (2.14). By its 
Green function (or fundamental solution) we mean the (generalized) function 
G e a Jx, t) which, being the formal solution of (2.2) corresponding to ip(x) = 
5(x) (the Dirac delta function), allows us to represent the solution of the 
Cauchy problem by the integral formula 



u. 



/+oo 
G e a ^,t)<p{x-i)di. (2.15) 
-oo 



It is straightforward to derive from (2.2) and (2.14) the composite Fourier- 
Laplace transform of the Green function by taking into account the Fourier 
transform for the Riesz-Feller space- fractional derivative, see (2.4), and the 
Laplace transform for the Caputo time- fractional derivative, see (2.12). We 
have (in an obvious notation) 

Gi^K, s) = f G*>, s) - S- 1 , (2.16) 

so that 

<p^ = sP+^m ' s>0 ' KGR (2 ' 17) 

In the special case 6 = we get 

^ (M) = TOF ' s>0 ' kgr - (2 - 18) 



Liouville fractional integral, is here denoted as tD° f(t) . We have, see e.g. [52] . 

/(r)dr' 



r(i - p) J (t - t) 



0</3< 1. 



It turns out that 

t D2 f{t) = t D? [f(t) - /(0+)] = t D ri f{t) - f(0+) Y^J) > < P < 1 • 

The Caputo fractional derivative, practically ignored in the mathematical treatises, 
represents a sort of regularization in the time origin for the Riemann-Liouville fractional 
derivative and satisfies the relevant property of being zero when applied to a constant. 
For more details on this fractional derivative (and its extension to higher orders) we refer 
the interested reader to Gorenflo and Mainardi [17] . 
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By using the known scaling rules for the Fourier and Laplace transforms, 
we infer directly from (2.17) (without inverting the two transforms) the 
following noteworthy scaling property of the Green function, 

G%(x,t) = t~^K e a ^(x/t^ a ) . (2.19) 

Here x/t^' a acts as the similarity variable and as the reduced Green 

function. Using (2.19) and the initial condition G e a o(x, + ) = 5{x) , we note 
that 

/+oo r+oo 
G e a ^(x,t)dx = / K e a>p {x) dx = 1 . (2.20) 
-oo J —oo 

In the case of the standard diffusion equation (2.1) the Green function is 
nothing but the Gaussian probability density with variance a 2 = 2t , namely 

G° 21 (x,t) = -± r t-^e- x2 /( 4t l (2.21) 

In the general case, following the arguments by Mainardi, Luchko & Pagnini 
[38] . we can prove that G^ t/3 (x,t) is still a probability density evolving in 
time. In the next section we summarise some results from [38) . 



3 The Green function for space-time fractional 
diffusion 

For the analytical and computational determination of the reduced Green 
function, from now on we restrict our attention to x > because of the 
symmetry relation ^(—x) = K~ 9 p(x) . Mainardi, Luchko & Pagnini [38] 
have provided (for x > 0) the Mellin-Barnes integral representation 



1 ! p+ioc r (±) r (i_i) r (i_ a ) 



K%(x) = / —x b ds, 

ax2mJ 1 - ioo r(l - §s)r(ps)r(l -pa) 



(3.3) 



a 



2a 

where < 7 < min{a,l}. Following [3H], we note that the Mellin- 
Barnea^l integral representation allows us to construct computationally the 



The names Mellin and Barnes refer to the two authors, who in the first 1910's 
developed the theory of these integrals using them for a complete integration of the 
hypergeometric differential equation. We note that, as pointed out in [TT] (Vol. 1, Ch. 
1, §1.19, p. 49), these integrals were first used by S. Pincherle in 1888. For a revisited 
analysis of the pioneering work of Pincherle we refer to [35] . 



7 



fundamental solutions of (2.2) for any triplet {a, 9, /?} by matching their 
convergent and asymptotic expansions. Readers acquainted with Fox H 
functions can recognize in (3.3) the representation of a certain function of 
this class, see e.g. [281 HH [52j 123 ES [61]. Unfortunately, as far as we know, 
computing routines for this general class of special functions are not yet 
available. 

Let us now point out the main characteristics of the peculiar cases of 
strictly space fractional diffusion and strictly time fractional diffusion, for 
which the non-negativity of the corresponding reduced Green functions is 
known. For (3 = 1 and < a < 2 (strictly space fractional diffusion) we 
have 

o , s ft, x 11 n +io ° T(s/a)T(l - s) „ , . , 

= L e a (x) = —— 1 ' ') } -x s ds, *>0, 3.6 

ax 2m J 7 _ioo 1 (p s) 1 (1 — ps) 

with < 7 < min{a, 1} , where L° a (x) denotes the class of the strictly stable 
(non-Gaussian) densities exhibiting heavy tails (with the algebraic decay 
oc |x|~( Q+1 )) and infinite variance. For a = 2 and < (3 < 1 (strictly time 
fractional diffusion) we have 

K°(x) = -M 0/2 (x) = / p n i \ x s ds, x>0, (3.7) 

Q ' U ' 2 W2V ' 2x2vrii 7 - i0 o T(l-(3s/2) ' v ' 

with < 7 < 1, where ^M^ 2 W denotes the class of the Wright- 
type densities exhibiting stretched exponential tails and therefore finite 
variance. The corresponding Green function evolves in time with the 
variance proportional to t@ . Mathematical details can be found in |38j ; 
for further reading we refer to Schneider [56] for stable densities, and to 
Gorenflo, Luchko & Mainardi [16] for the Wright-type densities. For the 
special case a = (3 < 1 referred in [38] as neutral diffusion, we obtain from 
(3.3) an elementary (non-negative) expression 

<«(*)= = — r(f)r(i-f) g 



ax 2m J-y-ioo T(ps)T(l — ps) 

(3.8) 

1 1 n+ iOC sin(Trps) s , 1 x^ 1 sin[^(a - 9)} 

x s ds = - _ - J\ , ^>0, 



ax2mjj-i QO sin(-7rs/a) it 1 + 2x a cosf^a — 6)] + x 2a ' 

with < 7 < a , where N®(x) denotes a peculiar class of densities exhibiting 
a power-law decay oc |cc| - ( a+1 ) , which contains the well known (stable) 
Cauchy density (recovered for a = 1 and 6 = 0). 
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For the generic case of strictly space-time diffusion (0 < a < 2, < 
(3 < 1), including neutral diffusion, we can prove the non negativity of the 
corresponding reduced Green function in virtue of the identity, see [38], 

K e p(x) = a Hf -1 M (Hi L e a (x/0 ^ , < /? < 1 , x >0 . (3.9) 

Then, as a consequence of the previous discussion, for the strictly space-time 
fractional diffusion we obtain a class of probability densities (symmetric or 
non-symmetric according to 9 = or 9 ^ 0) which exhibit heavy tails with 
an algebraic decay oc |x| . Thus they belong to the domain of attraction 

of the Levy stable densities of index a and can be referred to as fractional 
stable densities, according to a terminology proposed by Uchaikin [60]. In 
Fig. 1 we exhibit some plots of the probability densities provided by the 
reduced Green function for some "characteristic" values of the parameters 
a, /3, and 9. These plots, taken from [38], were drawn for the values of 
the independent variable x in the range |x| < 5. To give the reader a 
better impression about the behaviour of the tails the logarithmic scale was 
adopted. 

As for the stochastic processes governed by the above probability 
distributions we can expect the following. For the case of non-Gaussian 
stable densities we expect a special class of Markovian processes, called 
stable Levy motions, which exhibit infinite variance associated to the 
possibility of arbitrarily large jumps {Levy flights), whereas for the case 
of Wright-type densities we expect a class of stochastic non-Markovian 
processes, which exhibit a (finite) variance consistent with slow anomalous 
diffusion. Finally, for the case of fractional stable densities, the related 
stochastic processes are expected to possess the characteristics of the 
previous two classes. Indeed, they are non-Markovian (being (3 < 1) and 
exhibit infinite variance associated to the possibility of arbitrarily large 
jumps (being a < 2). A way to realize (understand) all the above stochastic 
processes is to show sample paths and histograms of related random walk 
models. For random walks discrete both in space and in time we refer to 
our papers [EH EH E2 [23] . 



4 From CTRW to fractional diffusion 

Here we show how the space-time fractional diffusion equation (2.2) can 
be obtained from the master equation for a continuous time random walk 
or, equivalently, from the master equation describing a cumulative renewal 
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-5 -4 -3 -2 -1 1 2 3 4 S -5 -4 -3 -2 -1 1 2 3 4 5 



Figure 1: Probability densities (reduced Green functions) for some values of 
the triplet {a, 9, (3} 
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process, through an appropriate limit. As a matter of fact this limit will be 
carried out in the Fourier-Laplace domain, so the corresponding convergence 
is to be intended in a weak form, that is sufficient for our purposes. For 
the basic principles of continuous time random walk (simply referred to as 
CTRW), that was formerly introduced in Statistical Mechanics by Montroll 
and Weiss [36], see e.g. [301 H5l Wl\ [63] . of renewal processes, see e.g. fTOl [TBI 

E5H5ZJ. 

The CTRW arises by a sequence of independently identically distributed 
(iid) positive random waiting times Ti, T2, . . . , each having a pdf ip(t) , t > 
0, and a sequence of iid random jumps Xi,X2,X%, .... in R, each having 
a pdf w(x) , x G R. Setting t = , t n = T± + T% + ■ ■ ■ T n for n € N , 
< t\ < t2 < ■ ■ ■ , the wandering particle starts at point x = in instant 
t = and makes a jump of length X n in instant t n , so that its position is 



x = for < t < Ti = ti , 

x = S n = X\ + X 2 + . . . X n , for t n <t < t n+1 . 

An essential assumption is that the waiting time distribution and the jump 
width distribution are independent of each other. It is well known that this 
stochastic process is Markovian if and only if the waiting time pdf is of the 
form ip(t) = mexp(-mt) with some positive constant m {compound Poisson 
process), see e.g. [13]. Then, by natural probabilistic arguments we arrive at 
the master equation for the spatial pdf p(x,t) of the particle being in point 
x at instant t , see [251 SOI EH] , 



p{x,t) = 6(x)V(t)+ / i>(t-t') 
Jo 



w(x — x) p(x , t) dx' 



dt' , (4.1) 



in which 5(x) denotes the Dirac generalized function, and, for abbreviation, 
^f(t) = J* t °° ip(t') dt' , is the probability that at instant t the particle is still 
sitting in its starting position x = 0. For this reason the function *&(t) 
is usually referred to as the survival probability; in the Markovian case it 
reduces to the exponential function *&(t) = exp(— mt) . Actually, p(x,t) as 
containing a point measure, is a generalized pdf, but for ease of language 
we omit the qualification "generalized". Clearly, (4.1) satisfies the initial 
condition p(x, + ) = 6(x) . 

It is customary (and convenient for our purposes) to consider the master 
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equation (4.1) in the Fourier-Laplace domain^ where it reads 

p(n, s) = $(s) + ip(s) w(k) p(k, s) , (4.2) 

whence, 

p(K,s) = — = ^ = — . (4.3) 

1 — w(k) ip(s) s 1 — w(k) V'(s) 

We will henceforth assume that in our continuous time random walk the 
jump width pdf w(x) is an even function (w(x) = w(—x)) and has a finite 
second moment (variance) or exhibits the asymptotic behaviour w(x) ~ 
^ | x |-(a+i) with some a , < a < 2 , for \x\ — > oo , and the waiting time pdf 
vp(t) has a finite first moment (mean) or exhibits the asymptotic behaviour 
ip(t) ~ ct~^ +1 ) with some (5 , < /3 < 1 , for t — > cxd , where 6 and c are 
positive constants. 

Our aim is to derive from the master equation (4.1), by properly rescaling 
the waiting times and the jump widths and passing to the diffusion limit, the 
space-time fractional diffusion equation. By our derivation of (2.2) from (4.1) 
we de-mystify the often asked- for meaning of the time fractional derivative in 
the fractional diffusion equation. In plain words, the fractional derivatives in 
time as well as in space are caused by asymptotic power laws and well-scaled 
passage to the diffusion limit. 

Scaling is achieved by making smaller all waiting times by a positive 
factor r , and all jumps by a positive factor h . So we get the jump instants 

t n (T) = + tT 2 + . . . + rT n for n G N , (4.4) 

and the jump sums, 

S (h) = , S n {h) = hXi + hX 2 + . . . + hX n for n G N . (4.5) 

The reduced waiting times rT n all have the pdf ip T (t) = ip(t/r)/T, t > 
0, and analogously the reduced jumps hX n all have the pdf Wh(x) = 
w(x/h)/h , iGR. Readily we see 

Tpr(s) = Ip(sT) , W h (K) = w(kK) . (4.6) 



9 It was in this domain that originally in 1965 Montroll and Weiss [46] derived their 
celebrated equation for the CTRW in Statistical Mechanics. However such equation can 
be derived by simply considering a random walk subordinated to a time renewal process 
as noted by us in |24| and by Baeumer and Meerschaert in [3]. 
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Replacing in (4.1) if)(t) by ip T (t) , V(t) by * T (i) = J™ ip T (t') dt' , w(x) by 
Wh(x) , p(x,t) by ph jT (x,t) we obtain the rescaled master equation which in 
the Fourier-Laplace domain reads as 

P/i,r(K, S) = 1 h ^r(s) t«/i(«) PhA K l S ) ' ( 4 - 7 ) 

s 

whose solution is 

S 1 - W h (K) 1p T {s) 

To proceed further we assume the probability densities w(x) and ip(t) of 
the jumps X n and the waiting times T n to meet the asymptotic conditions 
of the following Lemma 1 and Lemma 2, respectively, herewith recalled from 
where the interested reader can find the proofs. 

The first Lemma is a modified specialisation of Gnedenko's theorem in 
see also [8]. It was already used by us, but not formally called a Lemma, 
in [20]. The second Lemma can be obtained by aid of a corollary in Widder's 
book 1641. 



Lemma 1 



Assume w{x) > , w{x) = w(—x) for x S R, w(x) dx = 1 , and either 

r+oo 

a 2 := / x 2 w(x) dx < oo (4.9) 



oc 



(relevant in the case a = 2) or, with b > and some a £ (0, 2) , 

w {x) = (b + e(\x\))\x\-( a+1 K (4.10) 

In (4.10) assume e(|x|) bounded and O (\x\~ v ) with some r/ > as \x\ — > oo . 
Then, with a positive scaling parameter h and a scaling constant 

f V > if a = 2 , 

{ r(a + 1) sin(avr/2) ' if °<«< 2 ' 
we have, for each fixed k € R , the asymptotic relation 

w{ K h) = 1- fi(\ K \h) a +o{h a ) for fc->0. (4.12) 
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We note that (4.12) holds trivially if n = since w(0) = 1 
Lemma 2 



Assume tp(t) > for t > , J °° ifi(t) dt = 1 , and either 

p := / tif)(t)dt < oo (4.13) 
j o 

(relevant in the case (3=1), or, with c > and some (3 G (0, 1) , 

V>(i) ~ ct~ (/3+1) for t -> oo . (4.14) 
Then, with a positive scaling parameter r and a scaling constant 

fp, if /9 = 1, 

A = I cr(l-/3) ( 4 - 15 ) 
[ U W , if 0</3<l, 

we have, for each fixed s > , the asymptotic relation 

^( sr ) = 1 - A (srf + o(t P ) for r^O. (4.16) 
We note that (4.16) holds trivially if s = since tjj{0) = 1. 

Eq. (4.8) then becomes asymptotically 



By imposing the scaling relation 



Ph,A K > s ) ~ x g a , 7 q, I , a , fo r k,r->0. (4.17) 



\T f3 = fih a , (4.18) 

the asymptotics (4.17) yields 

P/t , T ("^)^ si9 + Na - ( 4 -f9) 

Hence, in view of (2.18), 

^(K,a)->G^(/c,«) for /i,r^0, (4.20) 

under condition (4.18). Then, the asymptotic equivalence in the space- 
time domain between the master equation (4.1) after rescaling and the 
fractional diffusion equation (2.2) with 8 = and the initial condition 
u(x,0 + ) = 5(x) is provided by the continuity theorem for sequences of 
characteristic functions after having applied the analogous theorem for 
sequences of Laplace transforms, see e.g. [13]. Therefore we have convergence 
in law or weak convergence for the corresponding probability distributions. 
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5 Simulations 



By aid of the results of Section 4 we can produce approximate particle paths 
for space-time fractional diffusion in the spatially symmetric case 8 = 
of (2.2). To this end, we require, for given a and (3, a jump width pdf 
w(x), obeying Lemma 1, and a waiting time pdf ip(t), obeying Lemma 2. 
Natural choices are the corresponding symmetric stable density of index 
a, i.e.w(x) = L^(x) (0 < a < 2) and, following [3D], the corresponding 
Mittag-Lefher type function 

^(t) = -jE p (-tP), 0</3< 1, (5.1) 

where 

W-Efj^- f>0, ,€C, (5.2) 

denotes the (entire) transcendental function, known as the Mittag-Lemer 
function of order (3 [UJ (Vol. 3, Ch 18, pp. 206-227). This function, which 
is a natural generalization of the exponential to which it reduces as (3 = 1 , 
is playing a fundamental role in the applications of fractional calculus, see 
e.g.[T71 [37] ■ As has been shown in [3D] , see also [28] [29], this choice of 
waiting-time density leads from the master equation (4.1) to the equation 

r+oo 

tD* p(x,t) = —p(x,t) + w(x — x')p(x',t)dx', p(x, + ) = 5(x) , (5.3) 

J — oo 

from which, by an appropriately scaled limiting process (analogous to that 
of Section 4), the fractional diffusion equation (2.2) with u(x,0 + ) = 5(x) , 
can be deduced, see [251. Observe that (5.3) in the particular case (3 = 
1 reduces to the well-known Feller-Kolmogorov equation for a compound 
Poisson process, in accordance with Ei(—t) = exp(— t) . 

Still some work must be invested in the inversions of the cumulative 
function W(x) = J^ 00 w(x') dx' and the survival probability ^(t) = 
J t °° ip{t') dt' , which here is 

$(t) = Epi-t 13 ) , t>0, 0</3<l. (5.4) 

In Fig. 2 we exhibit plots of *&(t) versus time for some values of j3 E (0, 1] 
from which we can get insight into the different behaviour for < [3 < 1 
(fast decay for short times and slow decay for long times) and for (3 = 1 
(exponential decay). 
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Figure 2: The survival probability V(t) = E^-t 13 ) for (3 = 0.25, 0.50, 0.75, 1 

These inversions are required by the standard Monte-Carlo procedure 
of generating the corresponding jump-widths and waiting-times from [0, 1] 
- uniformly distributed (pseudo-) random numbers. A. Vivoli in his thesis 
|62j has described in detail how all this can be done and has carried out 
several case studies of which we show (here) three samples for CTRW's, just 
to convey a visual impression on the structure of such processes, see Fig. 3. 
In these samples we have a = 2 so the jump density is a Gaussian, whereas 
(3 = 1,0.75,0.50. 

Observe in Fig. 3 the striking contrast between the first graph and the 
other two. In the case [3 = 1 we have ^(t) = exp(— t) which results in long 
waiting times occurring rarely (the mean waiting time being finite!). So, we 
get a good approximation of Brownian motion, (2.2) reducing to (2.1). For 
< (3 < 1, however, the Mittag-Leffler function exhibits a power-law decay, 
namely 

T , x ^ / Ax sin(/?7r) T(3) . . 

V(t) = Ef,(-tP) — t^oo. (5.5) 

As a consequence, we have a distinctly visible preponderance of long waiting 
times (the mean waiting time being infinite!). 

As our emphasis in this paper is on waiting times (relevant in CTRW's) 
we should say that the essential aspect is the asymptotic behaviour for 
t — > oo of the corresponding probability densities, namely, according to 
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Figure 3: Sample paths for CTRW's with a = 2, 9 = and = 1, 0.75, 0.50 
(from top to bottom) 



Lemma 2, their decay like ci~^ +1 ) (0 < (3 < 1) which implies for the survival 
probability a decay like (c//5)t _/3 . This is true, of course, for the Mittag- 
Leffler waiting-time distributions used here, see (5.1) and (5.5). However, in 
the interest of easy inversion of ^f(t), it is advantageous to look for simpler 
suitable functions. One such function (which is more easily invertible) is 

** (t) = i + r(i-flf/> ' t - ' 0</3<1 ' ( 5 - 6 ) 

so that 

**(*)(*) r la-) t^co. (5.6) 

7T 

Happily this function shares with the function ^(t) = Ea(—1r) the desirable 
property of complete monotonicity in t > 0. Furthermore we note that 
the functions ^(t) and VP*(i) share the same order of asymptotics for t —* + 
(albeit with a different coefficient). In fact we find as t — > + , 



In a forthcoming paper [26] we will describe in more detail our methods of 
simulation and investigate their quality. In the interest of long-time (or, 
because of self-similarity, near-the-limit) simulations, it is highly desirable 
that such fast methods are developed. 



6 Conclusions 



In this paper we have surveyed the general theory of the one-dimensional 
space-time fractional diffusion equation and have presented representation 
of its fundamental solutions (the probability densities) in terms of Mellin- 
Barnes integrals. Then, we have outlined how, in the spatially symmetric 
case, this equation can be obtained by a limiting process from a master 
equation for a continuous time random walk via properly scaled compression 
of waiting times and jump widths. For the strictly space and/or time 
fractional cases ({0 < a < 2, 0< (3 < 1}), it suffices to assume asymptotic 

d n 

10 Complete monotonicity of a function fit), t > 0, means (— l) n - — f(t) > 0, 

(n = 0,1,2,...), a characteristic property of exp(— t). For the Bernstein theorem this 
is equivalent to the representability of f(t) as (real) Laplace transform of a given non- 
negative (ordinary or generalized) function. For more information, see e.g. [4] (pp. 61-72), 
P3] (pp. 335-338), [31] (pp. 162-164) and [44]. 
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power laws 6|x| _(a+1 ) as \x\ — > oo for the jump width pdf and ct^^ +1 ^ as 
t — ► oo for the waiting time pdf. For the compression factors /i in space and 
t in time we require a scaling relation of the kind A h a where X,fi 

are given positive constants. Here we have limited ourselves to show sample 
paths for some cases of the time fractional diffusion processes (the jump 
width pdf is Gaussian) , referring for more comprehensive numerical studies 
to a forthcoming paper. The theory can be generalized to more than one 
space dimension and to non-symmetric jump pdf's, likewise to probability 
distribution functions instead of densities for the jump widths and waiting 
times, but, in order to avoid too cumbersome notations and calculations, let 
us just hint here to such possibilities. 
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